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Abstract 

By performing axisymmetric liydrodynamic simulations of core-collapse supernovae with spectral neu- 
trino transport based on the isotropic diffusion source approximation scheme, we support the assumption 
that the neutrino-heating mechanism aided by the standing accretion shock instability and convection can 
initiate an explosion of a 13 Mq star. Our results show that bipolar explosions arc more likely to be associ- 
ated with models which include rotation. We point out that models, which form a north-south symmetric 
bipolar explosion, can lead to larger explosion energies than for the corresponding unipolar explosions. 
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1. Introduction 

Core-collapse supernovae have long attracted the at- 
tention of astrophysicists because they have many facets 
playing important roles in astrophysics. They herald the 
birth of neutron stars and black holes; they are a major 
site for nucleosynthesis; they influence galactic dynamics; 
they trigger further star formation and they are prodigious 
emitters of neutrinos and gravitational waves. Despite 
rigorous theoretical studies for more than 40 years, the 
details of the explosion mechanism have been obscured 
under the thick veils of massive stars. 

For more than two decades, the neutrino-heating mech- 
anism (Wilson 1985; Bcthe & Wilson 1985), relying on 
the energy deposition via neutrinos behind the stalled 
shock, has been supposed as the most promising sce- 
nario. However, one important lesson we have learned 
from the work of Liebendorfer et al. (2001); Rampp & 
Janka (2002); Thompson et al. (2003); Sumiyoshi et al. 
(2005) is that the neutrino heating, albeit with the best 
input physics and numerics to date, fails in spherical sym- 
metry (ID) (see, however, Kitaura et al. 2006). 

Pushed by supernova observations of the blast mor- 
phology (e.g., Wang et al. 2001; Tanaka et al. 2007), it 
is now almost certain that the breaking of the spheri- 
cal symmetry is the key to the supernova puzzle. The 
multi-dimensional (multi-D) hydrodynamic motion asso- 
ciated with convectivc overturn in the postshock region 
(Herant et al. 1994; Burrows et al. 1995; Janka & Mueller 
1996; Fryer & Warren 2002; Fryer & Warren 2004) and 
the recently identified standing accretion shock insta- 



bility (SASI) (e.g., Blondin et al. 2003; Scheck et al. 
2004; Ohnishi et al. 2006; Foglizzo et al. 2007; Murphy 
& Burrows 2008; Iwakami et al. 2008; Iwakami et al. 
2009; Guilet ct al. 2009), arc expected to help the 
neutrino-driven explosion mechanism. This is because the 
sojourn time of the accreting material in the gain region 
can be longer than in the ID case, which enhances the ef- 
ficiency of the energy deposition behind the stalled shock. 

In fact, several explosion models have been reported 
recently in simulations that include multi-D effects that 
increase the neutrino heating (Buras et al. 2006; Marek & 
Janka 2009; Bruenn et al. 2009). Based on the long-term 
two-dimensional (2D) simulations with one of the best 
available neutrino transport approximations, Buras et al. 
(2006) firstly report an explosion for the non-rotating low- 
mass (11.2 Mq) progenitor of Woosley et al. (2002). Due 
to the compactness of the iron core 1.26Mq) with its 
steep outer density gradient, the explosion is initiated at 
^ 300 ms after core bounce. This is much earlier than 
in Marek & Janka (2009), who observe the delayed on- 
set of the explosion ~ 600 ms for a 15M0 progenitor of 
Woosley & Weaver (1995) with a moderately rapid ro- 
tation imposed. Although the explosion mechanism by 
neutrino-heating is very plausible, there are other pos- 
sible mechanisms, in which the magnetohydrodynamic 
(MHD) mechanism is included (see references in Kotake 
et al. 2006; Obcrgaulinger et al. 2006; Burrows et al. 
2007b; Takiwaki et al. 2009). Another suggested mech- 
anism relics on acoustic energy deposition via oscillating 
protoneutron stars (PNSs), which has been discovered by 
a series of 2D multi-energy flux-limited-diffusion transport 
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simulations (Burrows et al. 2006; Burrows et al. 2007a). 
Although the additional energy input from sound appears 
to be robust enough to explode even the most massive 
progenitors (Burrows et al. 2007a), it remains a matter of 
vivid debate and has yet to be confirmed by other groups. 
Also exotic physics in the core of the PNS may have a 
potential to trigger explosions (e.g., Sagert et al. 2009). 

In this Letter, we present axisymmetric explosion mod- 
els for a 13 Mq progenitor model of Nomoto & Hashimoto 
(1988) in support of the theory that neutrino-heating 
aided by multi-D effects is able to cause supernova ex- 
plosions. We choose the progenitor with a smaller iron 
core (~ I.2OM0), anticipating an explosion since the pro- 
genitor mass lies between 11.2 Mq (Buras et al. 2006) 
and 15Mq (Marek & Janka 2009). We perform 2D 
core-collapse simulation with spectral neutrino transport 
by the isotropic diffusion source approximation (IDSA) 
scheme currently developed by Liebendorfer et al. (2009). 
By comparing four exploding models with and without 
rapid rotation to one non-exploding ID model, we point 
out that models that produce a north-south symmetric 
bipolar explosion can lead to larger explosion energies 
than for the corresponding unipolar explosions. Our re- 
sults show that the explosion geometry is more likely to 
be bipolar in models that include rotation. 

2. Numerical Methods and Models 

Our 2D simulations are performed using a newly de- 
veloped code which implements spectral neutrino trans- 
port using the IDSA scheme (Liebendorfer et al. 2009) 
in a ZEUS-2D code (Stone & Norman 1992). Following 
the spirit of the so-called ray-by-ray approach, the IDSA 
scheme further splits the neutrino distribution into two 
components, both of which are solved using separate nu- 
merical techniques. Although it does not yet include 
heavy lepton neutrinos such as v^^Vt {vinVr) and the in- 
elastic neutrino scattering with electron, the innovative 
approach taken in the scheme saves a significant amount of 
computational time compared to the canonical Boltzmann 
solvers (see Liebendorfer et al. 2009 for more details). 
Expecting a bigger chance to produce explosions (Marek 
& Janka 2009), we employ the soft equation of state (EOS) 
by Lattimer & Swesty (1991) with a compressibility mod- 
ulus oi K = 180 MeV. The self gravity is implemented by 
solving the Poisson equation by the Modified Incomplete 
Cholesky Conjugate Gradient (MICCG) method (Kotake 
et al. 2003), but without relativistic corrections. 

The simulations are performed on a grid of 300 loga- 
rithmically spaced radial zones up to 5000 km. To test 
the sensitivity with respect to angular resolution, the grid 
is varied to consist of 64 or 128 equidistant angular zones 
covering < < tt. For the neutrino transport, we use 20 
logarithmically spaced energy bins reaching from 3 to 300 
MeV. 

All supernova calculations in this work are based on 
the 13Mq model by Nomoto & Hashimoto (1988). The 
computed models are listed in the first column of Table 
1, in which one calculation (model M13-1D) is conducted 
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Fig. 1. Time evolution of Models M13-1D and M13-2D, vi- 
sualized by mass shell trajectories in thin gray and orange 
lines, respectively. Thick lines in red (for model M13-2D) 
and black (model M13-1D) show the position of shock waves, 
noting for 2D that the maximum (top) and average (bottom) 
shock position are shown. The red dashed line represents the 
position of the gain radius, which is similar to the ID case 
(not shown). 

in spherical symmetry. Other models are 2D simulations 
with or without rotation (indicated by rot) with different 
numerical resolution in the lateral direction (64 or 128, 
denoted by "hr" (high resolution) in Table 1). For the ro- 
tating models, we impose rotation on the progenitor core 
with initially a constant angular frequency of fJo = 2 rad/s 
inside the iron core with a dipolar cut off (oc r^^) outside, 
which corresponds to /? ~ 0.18% with /? being the ratio of 
the rotational to the gravitational energy. This rotation 
rate is fairly large and may lead to a spiral mode of the 
SASI (Yamasaki & Foglizzo 2008). In addition, this strong 
rotation may induce a strong magnetic field due to wind- 
ing and the magneto-rotational instability and produce a 
jet-like outflow (Kotake et al. 2006). Although these ef- 
fects could modify the dynamics of the postbounce phase, 
the approximate treatment in this study (axisymmetry 
without magnetic fields) does not allow us to investigate 
them in this article. 

3. Results 

Figure 1 depicts the difference between the time evolu- 
tions of model M13-1D (thin gray lines) and model M13- 
2D (thin orange lines), visualized by mass shell trajecto- 
ries. Until ~ 100 ms after bounce, the shock position of 
the 2D model (thick red line) is similar to the ID model 
(thick black line). Later on, however, the shock for model 
M13-2D docs not recede as for M13-1D, but gradually ex- 
pands and reaches 1000 km at about 470 ms after bounce. 
Comparing the position of the gain radius (red dashed 
line) to the shock position of M13-1D (thick black line) 
and M13-2D (thick red line), one can see that the ad- 
vcction time of the accreting material in the gain region 
can be longer in 2D than ID. This longer exposure of 
cool matter in the heating region to the irradiation of 
hot outstreaming neutrinos from the PNS is essential for 
the increased efficiency of the neutrino heating in multi-D 
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Fig. 2. Snapshot of the distribution of entropy (left half) and the ratio of the advcction to the heating timescale (right half) for 
models of M13-1D (left) and M13-2D (right) at 200 ms after bounce. 



models. 

A more detailed analysis of the timescale is shown in 
Figure 2. The right-half shows Tadv/'Hicat, which is the 
ratio of the advection to the neutrino heating timescale. 
For the 2D model (right panel), it can be shown that 
the condition of Tadv/'Ticat ^ 1 is satisfied behind the as- 
pherical shock, which is deformed predominantly by the 
SASI, while the ratio is shown to be smaller than unity in 
the whole region behind the spherical standing accretion 
shock (left panel: ID). Note that Thcat is estimated locally 
by ebind/Qi/, where Cbind is the local specific binding en- 
ergy (the sum of internal plus kinetic plus gravitational 
energies) and Qi, is the specific heating rate by neutrinos, 
and that Tadv is given by [r~rgain{0)]/\vr{r,e)\, where rgai„ 
is the gain radius and Vr is the radial velocity. Comparing 
the left-half of each panel, the entropy for the 2D model is 
shown to be larger than for the ID model. This is also the 
evidence that the neutrino heating works more efficiently 
in multi-D. 

We now move on to discuss models with rotation. Both, 
for model M13-rot and its high resolution counterpart, 
model M13-rot-hr, we obtain neutrino-driven explosions 
(see, tiooo and .Bdia in Table 1). The rapid rotation chosen 
for this study mainly affects the explosion dynamics in the 
postbounce phase, which we will discuss in the following. 

For the rotating model, the dominant mode of the shock 
deformation after bounce is almost always the £ = 2 mode 
although the £=1 mode can be as large as the £ — 2 mode 
when the SASI enters the non-linear regime (^ 200 ms 
after bounce) . In contrast to this rotation- induced £ = 2 
deformation, the £ = 1 mode tends to be larger than the 
£ = 2 mode for the 2D models without rotation in the satu- 
ration phase. As shown in Figure 3, this leads to different 
features in the shock geometry, namely the preponder- 
ance of the unipolar explosion for the 2D models without 
rotation (left), and the bipolar (north-south symmetric) 
explosion with rotation (right). 

Since it is impossible to calculate precise explosion ener- 
gies at this early stage, we define a diagnostic energy that 
refers to the integral of the energy over all zones that have 



a positive sum of the specific internal, kinetic and grav- 
itational energy. Figure 4 shows the comparison of the 
diagnostic energies for the 2D models with and without 
rotation. Although the diagnostic energies depend on the 
numerical resolutions quantitatively, they show a contin- 
uous increase for the rotating models. The diagnostic en- 
ergies for the models without rotation, on the other hand, 
peak around 180 ms when the neutrino-driven explosion 
sets in (see also Figure 1), and show a decrease later on. 
With values of order 10*^ erg it is not yet clear whether 
these models will also eventually lead to an explosion. 

The reason for the greater explosion energy for models 
with rotation is due to the bigger mass of the explod- 
ing material. This is because the north-south symmetric 
{£ = 2) explosion can expel more material than for the 
unipolar explosion. In fact, the mass enclosed inside the 
gain radius is shown to be larger for the rotating models 
(e.g.. Table 1). The explosion energies when we termi- 
nated the simulation are less than ^ 10^° erg for all the 
models. For the rotating models, we are tempted to spec- 
ulate that they could become as high as ^ 10^^ erg within 
the next 500 ms by a linear extrapolation. However, in 
order to unquestionably identify the robust feature of an 
explosion in the models, a longer-term simulation with 
improved input physics would be needed. 

Our numerical results are qualitatively consistent with 
the results of Marek & Janka (2009) in the sense that in a 
relatively early postbounce phase the model with rotation 
shows a more clear trend of explosion than for the non- 
rotating models. 

4. Summary and Discussion 

By performing 2D core-collapse simulations of a 13 Mq 
star with spectral neutrino transport via the isotropic dif- 
fusion source approximation, we have found a strong de- 
pendence of the expansion of the shock radius and the 
likelihood for an explosion on the initial rotation rate. In 
all cases the shock is driven outward by the neutrino heat- 
ing mechanism aided by multi-D effects such as the SASI 
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Fig. 3. Snapshots of the density (left half) and the entropy (right half) for models M13-2D (left) and M13-rot (right) at the epoch 
when the shock reaches 1000 km, corresponding to ~ 470 ms after bounce in both cases. 

be affected so much because they are produced mainly by 
the hydrodynamic interplay of the SASI and the rotation. 
The ray-by-ray approach may lead to the ovcrcstimation 
of the directional dependence of the neutrino anisotropics 
(see discussions in Marek & Janka 2009). On the other 
hand, the lateral neutrino emission and the enhanced 
heating near the polar regions, such as from the oblately 
deformed protoneutron star due to rapid rotation (e.g., 
Kotake et al. 2003), could be underestimated. Apparently 
the full-angle transport will give us the correct answer 
(e.g., Ott et al. 2008). In addition, due to the coordinate 
symmetry axis, the SASI develops preferentially along the 
axis, thus it could provide a more favorable condition for 
the explosion. As several exploratory simulations have 
been done recently (Iwakami et al. 2008; Scheidegger et al. 
2008; Iwakami et al. 2009), 3D supernova models are in- 
deed necessary. 

Bearing these caveats in mind, the role of rotation act- 
ing on the neutrino-driven explosions, is qualitatively new. 
Yet there remain a number of issues to be studied. We 
have to clarify the progenitor dependence and also in- 
vestigate the effects of rotation more systematically by 
changing its strength in a parametric manner (possibly 
with magnetic fields). It will be interesting to study the 
neutrino and gravitational-wave signals. This paper is a 
prelude for the forthcoming work that will clarify these 
issues one by one. 
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Fig. 4. Time evolution of the diagnostic energy versus post- 
bounce time for 2D models with and without rotation. 

and convection. We have shown the preponderance of an 
bipolar explosion for 2D models with rotation. We have 
pointed out that the explosion energy can become larger 
for models with bipolar explosions. 

The conclusion with respect to the effects of rotation 
obtained in this study differs from that of Marek & Janka 
(2009), who suggested that the rotation has a negative 
impact on the explosion. They obtained the expansion of 
the shock wave only for the rotating model (M15LS-rot). 
while the nonrotating model did not show an expansion 
due to the short simulation time (see Fig. 6 in their pa- 
per). Therefore they could not compare the expanding 
shock evolution in both the rotating and the nonrotat- 
ing cases so that their discussion is limited to the shock 
oscillation phase. 

Here it should be noted that the simulations in this 
paper are only a very first step towards more realistic 
supernova models (e.g., Marek & Janka 2009; Burrows 
et al. 2007a; Bruenn et al. 2009). The approximations 
adopted in this paper should be improved, for example 
the omission of heavy lepton neutrinos, the inelastic neu- 
trino scattering, and the ray-by-ray approach. The for- 
mer two, may act to suppress the explosion. However we 
think that qualitative effects induced by rotation will not 
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Models 



Dimension 



[rad/s] 



^1000 

[msl 



Ed 



la -'^-'gain 

erg] [Mq] 



M13-1D 


ID 




1 








M13-2D 


2D 




64 


470 


0.26 


0.017 


M13-rot 


2D & rotation 


2 


64 


480 


0.95 


0.067 


M13-2D-lir 


2D 




128 


420 


0.40 


0.018 


M13-rot-hr 


2D & rotation 


2 


128 


520 


0.78 


0.060 



SZq is the precollapse angular velocity. Ng represents the lateral grid number covering 
< 6 < TT. "hr (high resolution)" indicates the runs for Ng = 128. iiooo represents the 
time (measured after bounce) when the average shock radius becomes 1000 km. -Edia is 
the diagnostic energy defined as the total energy (internal plus kinetic plus gravitational), 
integrated over all matter where the sum of the corresponding specific energies is positive. 
Mgain is the mass inside the gain layer. The latter quantities are given at 450 ms postbounce. 
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